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ABSTRACT 

We use ensembles of high-resolution CDM simulations to investigate the shape and 
amplitude of the two point correlation function of rich clusters. The standard scale- 
invariant CDM model with SI = 1 provides a poor description of the clustering mea- 
sured from the APM rich cluster redshift survey, which is better fitted by models with 
more power at large scales. The amplitudes of the rich cluster correlation functions 
measured from our models depend weakly on cluster richness. Analytic calculations of 
the clustering of peaks in a Gaussian density field overestimate the amplitude of the 
N-body cluster correlation functions, but reproduce qualitatively the weak trend with 
cluster richness. Our results suggest that the high amplitude measured for the corre- 
lation function of richness class R > 2 Abell clusters is either an artefact arising from 
incompleteness in the Abell catalogue, or an indication that the density perturbations 
in the early universe were very non-Gaussian. 

Key words: Galaxies : Clustering ; Large-scale structure of the Universe. 



Cr)l INTRODUCTION 

0^ . . 

Investigations of the spatial distribution of rich clusters pro- 
vided some of the earliest evidence for large scale irregulari- 
^4ies in the Universe (Kiang & Saslaw 1969, Yu & Peebles 
Q 1969, Bogart & Wagoner 1973, Hauser & Peebles 1973). 
^ More recently, various redshift surveys of clusters selected 
^from the Abell catalogue (Abell 1958) have been used to es- 
timate the rich cluster two-point correlation function ^cc{i') 
(e.g. Bahcall & Soneira 1983, Klypin & Kopylov 1983, Post- 
man, Huchra & Geller 1992, hereafter PHG). Two new clus- 
ter catalogues have been constructed from digitised scans of 
photographic plates, and redshift surveys have been carried 
out by Dalton et al. (1992) (the APM Cluster Catalogue), 
hereafter DEMS and Nichol et al. (1992) (the Edinburgh- 
Durham Cluster Catalogue, or EDCC). New redshift surveys 
of X-ray selected samples of rich clusters are in progress. 

It has been known since the work of Hauser and Peebles 
(1973) that rich clusters of galaxies are more strongly clus- 
tered than galaxies, but the amplitude of the cluster two 
point correlation function, and its dependence on cluster 
richness are controversial. Taken at face value, the three- 
dimensional clustering of rich clusters clusters is consistent 
with a scaling relation 
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U{r) ~ 0.2(r/rfe)" 



(1) 



where {cc(''o) = 1 and dc is the mean intercluster separa- 
tion related to the mean space density ric by dc - 
Bahcall & Burgett (1986), Bahcall & West (1992), and Bah- 
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call & Cen (1992) (hereafter BC92) argue that equation (1) 
applies for clusters over the range dc = 30 h~^Mpc * to 
dc X 90 /t"^Mpc. 

Redshift surveys based on the APM Cluster Catalogue 
and EDCC have provided accurate estimates of the corre- 
lation length To for clusters of richnesses intermediate to 
those of i? = and R = 1 Abell clusters. These surveys 
indicate ro ~ 13-16 /t~"'Mpc for clusters with dc ~ 35- 
46 /t~"'Mpc, consistent with equation (1). However, these 
surveys are too small to provide accurate measurements on 
the clustering of richer clusters. The extrapolation of equa- 
tion (1) to richer clusters relies on redshift surveys based 
on the the Abell catalogue (Abell 1958) and its southern 
counterpart (Abell, Corwin & Olowin 1989). For example, 
Bahcall and Soneira (1983) find ro w 25 /t"^Mpc for Abell 
R > 1 clusters {dc fa 55 /t~^Mpc), while Peacock & West 
(1992) find ro = 21.1 ± 1.3 h'^Mpc for i? > 1 Abell clus- 
ters and ro = 45 ± 5 /t"^Mpc for i? > 2 Abell clusters 
(dc w 84 fe-^Mpc). 

However, it is unclear whether these results for Abell 
clusters are correct. Various authors have presented evi- 
dence that the Abell catalogues, which were constructed by 
scanning photographic plates by eye, are affected by incom- 
pleteness on the plane of the sky which enhance the clus- 
tering amplitude measured in three dimensions (Sutherland 
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1988; Soltan 1988; Sutherland & Efstathiou 1991; Efstathiou 
et al. 1992, hereafter EDSM). A particularly strong case for 
patchy incompleteness in the Abell i? > catalogue has 
been presented by EDSM, who compare the machine based 
APM survey with the PHG redshift survey of Abell clus- 
ters. The clusters in these surveys have comparable space 
densities, but the redshift-space correlation function of the 
APM sample is isotropic on large scales with a low am- 
plitude (ro ~ 13 /t~^Mpc), while the correlation function 
for the Abell clusters is highly anisotropic and has a much 
higher amplitude. The effect of patchy incompleteness on 
samples of richer Abell clusters is more controversial. Suther- 
land (1988), EDSM and others have argued that the corre- 
lation functions of the i? > 1 Abell catalogue may have been 
significantly overestimated whereas some authors, e.g. Bah- 
call and West (1992) and Peacock & West (1992), suggest 
that any biases are small. 

Equation (1) is therefore approximately correct for clus- 
ters with dc ~ 35 /t~^Mpc and agrees with results from ma- 
chine measured surveys. However the extrapolation of equa- 
tion (1) to higher richnesses relies exclusively on empirical 
results from the Abell catalogue. Since there is evidence for 
patchy incompleteness in the Abell sample, the richness de- 
pendence implied by equation (1) must be viewed with some 
skepticism. In this paper we give high weight to the machine 
based surveys of DEMS and Nichol et al (1992) , since these 
are known to be free of large redshift-space anisotropics of 
the clustering pattern, and we give low weight to the re- 
sults from the Abell catalogue for the reasons outlined in 
the previous paragraph. We refer the reader to the papers 
by Dalton et al (1992) and Efstathiou et al (1992) for more 
detailed discussions of the observations. 

In this paper, we use N-body simulations to investigate 
the clustering of rich clusters from Gaussian initial condi- 
tions. Our aim is to determine the shapes of the cluster 
correlations and how these depend on the amplitude of the 
mass fluctuations and on cluster richness. A number of au- 
thors have investigated these problems theoretically using 
the statistical properties of Gaussian density fields (see e.g. 
Kaiser 1984; Bardeen et al. 1986, hereafter BBKS; Bardeen, 
Bond & Efstathiou 1987; Holtzman & Primack 1993; Mann, 
Heavens & Peacock 1993), but it is not clear how well these 
analyses compare with the full non-linear problem. White 
et al. (1987) and Bahcall & Cen (1992) have determined the 
two-point cluster correlations from N-body simulations of 
CDM-like universes, but our computations sample a larger 
volume of space than these studies, and have higher spatial 
resolution than the particle-mesh simulations of Bahcall & 
Cen. 

The layout of this paper is as follows. The numerical 
simulations and the cluster selection algorithm are described 
in Section 2. The cluster correlations are discussed in Sec- 
tion 3 and are compared to an analytic model based on 
the statistical properties of Gaussian density fields in Sec- 
tion 4. Section 5 compares the APM cluster redshift survey 
(DEMS) with simulated APM cluster catalogues generated 
from the N-body models. Our conclusions are summarized 
in Section 6. 



2 NUMERICAL SIMULATIONS AND 
CLUSTER SELECTION 

2.1 Numerical simulations 

We carried out 30 simulations of CDM-like universes in cubi- 
cal volumes of comoving box length Ib = 300 /t"^Mpc usmg 
a particle-particle-particle-mesh (P^M) code ( Efstathiou & 
Eastwood 1981, Efstathiou et al. 1985). Each simulation fol- 
lowed the non-linear gravitational evolution of 10^ particles 
in a 256'^ mesh with a spatial resolution of ~ 80/t~^kpc. The 
high resolution of our simulations ensures that we can follow 
the formation and evolution of mass concentrations that can 
be identified with rich clusters of galaxies (see Section 2.2). 
The initial power spectrum of our models is given by: 



(where j/ = 1.13, a = 6.4/r h'^Mpc, b = 3.0/T h'^Mpc, 
c = 1.7/r /t~^Mpc). This power spectrum applies for scale- 
invariant CDM models with low baryon densities Q_b ^ 
and r = Qo/i (Bond & Efstathiou 1984). Thus the 'standard' 
CDM model with Qo = 1 and A ~ 0.5 (Davis et al. 1985) 
has r Si 0.5. Ensemble A consists of 10 simulations of the 
standard CDM model with T = 0.5. There is considerable 
evidence that large-scale clustering of galaxies and clusters 
is better described by a power spectrum with T Si 0.2 (Efs- 
tathiou, Sutherland & Maddox 1990, Efstathiou 1993). We 
have therefore run two additional ensembles each consisting 
of 10 simulations of spatially flat models with T = 0.2, using 
the same phases as the models in ensemble A to generate the 
initial particle displacements. Ensemble B has a cosmologi- 
cal constant, A = A/(3_ffo^) = (1 ~ ^^o) = 0.8, and ensemble 
C has r = 0.2 and Qo = 1. 

The parameters of the models are given in Table 1. A, 
is the amplitude of the initial linear fluctuations relative to 
the white-noise level at the Nyquist frequency of the parti- 
cle grid used in setting up the initial conditions (Efstathiou 
et al. 1985). The column labelled (crs), gives the initial rms 
amplitude of the linear mass fluctuations in spheres of ra- 
dius 8 /t~^Mpc. Each model was evolved until the value of 
as was equal to 1.0 according to linear theory. The expan- 
sion factor required to reach this value of ag from the start 
of the simulations is given by afmai- 



Table 1. Parameters for the CDM-like models under investiga- 
tion 



Ensemble 


r 


h 


A 


A. 


((Tg), 


(1 final 


A 


0.5 


0.5 


0.0 


0.70 


0.236 


4.24 


B 


0.2 


1.0 


0.8 


0.64 


0.311 


4.50 


C 


0.2 


1.0 


0.0 


0.40 


0.194 


5.15 



Each simulation took about two days on a Decstation 
5000/240 with 168Mbytes of memory. With the large mesh 
used here, the calculations are mesh dominated; towards the 
end of the calculations, when the particle distributions are 
clustered most strongly, computing the short-range forces 
takes only about 25% of the time required to compute the 
mesh forces. 
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2.2 Cluster selection 

To select candidate centres for clusters we found groups of 
particles linked together by less than 0.1 times the mean 
interparticle separation using the percolation algorithm de- 
scribed by Davis et al. (1985). We then searched for clusters 
by counting particles within radius Tc of the percolation cen- 
tres. We calculated the centre of mass of each cluster and 
deleted any cluster within Tc of a cluster of a larger mass. 
We then repeated the cluster finding by counting particles 
within Tc of the centre of mass of each cluster, recomputing 
the centre of masses and deleting overlapping clusters. This 
step was repeated twice to determine the final list of clusters. 
This algorithm is the same as that used by White, Efstathiou 
& Frenk (1993). Having generated a cluster catalogue from 
a simulation for a specified value of Tc, we order the clusters 
by mass. By applying a lower mass limit, we can generate 
samples of different space densities Wc and mean intercluster 
distance dc = ric^^^ ■ 

This algorithm, which finds mass concentrations within 
spherical volumes in three-dimensions, is similar to the al- 
gorithms used to identify clusters from two-dimensional cat- 
alogues of galaxies. Abell's (1958) algorithm is designed to 
find clusters of galaxies within a metric radius of Ra = 
1.5 /t~"'Mpc. Clusters of galaxies in the APM survey are 
found using an Abell-like algorithm with a smaller metric 
radius of 0.75 /t~^Mpc, but since galaxies closer to the clus- 
ter centre are given higher weight than those in the outer 
parts, the effective cluster radius in the APM sample is 
slightly smaller. We have generated cluster catalogues us- 
ing Tc = 1.5 /t~"'Mpc and 0.5 /t~"'Mpc. As we will show in 
the next section, our results are almost independent of the 
value of Tc- 

The simple prescription that we have adopted avoids 
having to assign galaxies to mass points, which would intro- 
duce uncertain parameters and further assumptions since 
the galaxy formation process is so poorly understood. For a 
given choice of cluster radius Tc our procedure results in a 
unique list of clusters ordered by mass. Provided, therefore, 
that there is an approximate monotonic relation between the 
masses of rich clusters and their luminosity, our prescription 
should provide an accurate match to observed samples lim- 
ited by cluster richness, even though we make no specific 
model for assigning galaxies to the mass. This is an impor- 
tant simplification which we feel is reasonable and neces- 
sary to make the problem tractable. One can imagine situa- 
tions where our assumptions would be incorrect, for example 
if cluster luminosities were modulated by some correlated 
physical mechanism in analogy with models of 'cooperative 
galaxy formation' as discussed, for example by Bower et al. 
(1993). However, such mechanisms are so poorly understood 
that they have generally been ignored when comparing the- 
oretical models to observations. 

Figure 1 shows clusters in slices of depth 75 /t~"'Mpc 
at the final output times (erg = 1) from a model from 
ensemble A and one with identical initial random phases 
from ensemble B. The open circles show clusters identified 
within fc = 1.5 /t~^Mpc and the filled points show clus- 
ters identified within Tc = 0.5 /t~^Mpc. These catalogues 
have dc = 35 /t~^Mpc, comparable to the mean intercluster 
distance of 7J > 20 clusters in the APM survey (EDSM). 
The two pictures are visually extremely similar despite the 
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Figure 1. A 75 /i~^Mpc slice through simulations with iden- 
tical initial random phases picked from (a) ensemble A, and (b) 
ensemble B. Clusters are selected with Tc = 0.5 h~^Mpc (dots) 
and Tc = 1.5 h~^M.pc (Open circles), and with dc (mean inter- 
cluster separation) = 35 h~^Mpc. The symbol sizes are not to 
scale. 
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differences in tfie power-spectra. Tfiere is aiso a good corre- 
spondence between tfie ciusters found witfiin tfie two searcfi 
radii, ff aii ciusters fiad identicai density profiies, our ciuster 
cataiogues wouid be independent of Tc- fn fact, tliis is nearly 
tlie case and tlie clusters found with rc = 0.5 /t~^Mpc have 
approximately one third of the mass of clusters identified 
with Tc = 1.5 /t~"'Mpc, as expected for clusters with near 
isothermal profiles M{r) oc r (see also Figure 1 of White et 
al. (1993)). 



3 THE TWO-POINT CLUSTER 
CORRELATION FUNCTION 

For each cluster catalogue we compute the two-point corre- 
lation function using a direct estimator: 



where Np is the number of cluster pairs in the radial bin of 
volume dV centred at r, Wc is the mean space density of the 
cluster catalogue, V is the volume of the simulation. We use 
all clusters in a simulation taking advantage of the periodic 
boundary conditions. 

For each simulation we compute fee as a function of 
cluster richness and epoch. In the next subsection, we dis- 
cuss the evolution of the correlation functions. We then in- 
vestigate how fee depends on the cluster selection radius Tc 
and on cluster richness. 

3.1 Evolution of the two-point correlation 
function 

In Figure 2, we pick subsamples of clusters with a mean 
separation dc = 35 /t~^Mpc, close to that of APM Tl. > 
20 clusters (c/. Figure 1). We plot the correlation functions 
at various output times with the values of erg given in the 
Figure. 

The amplitude of the cluster correlations is insensitive 
to the amplitude of fluctuations in the density field. This 
is expected because these clusters are much more highly 
clustered than the mass fluctuations, even when the latter 
achieve values of erg ~ 1. The insensitivity of the clustering 
on (7g is observed in our simulations even for poor clusters 
with dc ~ 10 /t"^Mpc. There has been some controversy 
over the amplitude of the mass fluctuations in the universe 
at the present epoch (see for example Couchman & Carlberg 
(1992), who propose an Q = 1 CDM universe with erg = 1). 
However, a recent analysis of the abundances and masses of 
rich clusters of galaxies gives 

(7g = (0.57 ± 0.05) Q^"^^'' (4) 

(White, Efstathiou & Frenk 1993) and in the rest of this pa- 
per we use equation (4) to select appropriate output times 
in the simulations. Note that erg in equation (4) is the value 
extrapolated to the present day using linear pertubation the- 
ory, and is extremely insensitive to the shape of the initial 
fluctuation spectrum. As Figure 2 shows, fortunately we do 
not have to specify erg very precisely to make accurate pre- 
dictions of the cluster correlations. For our low density en- 
semble, a best fit to the COBE data requires that erg be 



slightly lower: ~ 1.0. The 10 models we plot results for there- 
fore have this value at the present day. We have also run 
4 additional low density models with erg = 1.34, and find 
the results to be indistinguishable from the less dynamically 
evolved models. 

In Figure 2 we have plotted the data points for the APM 
clusters computed by EDSM. We show points only out to 
50 /t"^Mpc since the results at larger radii are consistent 
with zero (eg. compare figure (10) of EDSM with figure (1) 
of DEMS). Note further that we plot Poisson error bars on 
the data points which we argue in Section (5) underestimate 
the true errors by as much as 50%. Although the observed 
correlations were computed in redshift space, we show in 
Section 5 that the distortions between redshift space and 
real space in the simulations are small. The F = 0.5 curves 
lie below the APM data points at separations > 8 /t~^Mpc. 
The two F = 0.2 models give much better agreement with 
the observations. 



3.2 Changing the radius of cluster selection 

Figure 3 shows the effect on each of our models of increas- 
ing the radius used to define clusters in the simulations from 
re = 0.5 fe-^Mpc to 1.5 fe-^Mpc. Since the correlation func- 
tions are so insensitive to the amplitude of the mass fluc- 
tuations (Figure 2), we plot results for one value of erg for 
each ensemble. The cluster catalogues were limited in mass 
to give an intercluster separation of dc = 35 /t~^Mpc. 

As noted in Section 2, there is a good correspondence 
between cluster catalogues with the same value of dc but 
different values of Tc (see Figure 1). Accordingly, the cor- 
relation functions for Tc = 0.5 and 1.5 /t~^Mpc shown 
in Figure 3 are very similar. The correlation functions for 
Tc = 0.5 /t~^Mpc are systematically lower than those for 
Tc = 1.5 /t~"'Mpc by about 10-20%, but these differences 
are much smaller than the observational errors on estimates 
of fee from real cluster catalogues. The independence of fee 
on Tc applies for all of the cluster catalogues we have gener- 
ated, i.e. for dc in the range 10-70 /t~^Mpc. 

As explained in Section 2, we chose Tc = 0.5 /t~"'Mpc 
to model the selection of clusters in the APM survey, and 
Tc = 1.5 /t~^Mpc to model the selection of Abell clusters. 
Figure 3 would lead us to expect that Abell and APM cluster 
catalogues with the same value of dc would have nearly the 
same amplitude of {cc We plot {cc for the APM Tl > 20 
clusters, as in Figure 2. The filled stars show results for the 
redshift survey of i? > Abell clusters (PHG) corrected for 
redshift-space anisotropics by EDSM. The Abell sample has 
dc = 38 /t~^Mpc, which is nearly the same as the value for 
the APM sample. The results for the two catalogues are very 
similar, and are well described by the F = 0.2 models. There 
is a marginal tendency for the Abell data points to lie higher 
than those for the APM survey, and this is in the direction 
suggested by the N-body results. 

As mentioned in the introduction, the clustering in the 
PHG sample agrees with that measured in the APM sam- 
ple only after the Abell catalogue has been corrected for 
incompleteness on the plane of the sky. The two-point cor- 
relation function for the uncorrected R > catalogue has 
about twice the amplitude of the APM correlation function 
(see PHG and EDSM). Clearly, in the class of models con- 
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Figure 2. Time variation of the cluster correlation functions for clusters with dc = 35 /i~^Mpc selected with a search radius of 
Tc = 0.5 /i~^Mpc. Figure 2(a) shows ensemble A, i.e. standard CDM with Qq = 1 and h = 0.5; (b) shows ensemble B in which F = 0.2, 
Qo = 0.2 and A = 0.8; (c) shows ensemble C in which F = 0.2 and Qq = 1- The error bars on the simulation curves show one standard 
deviation on the mean and were computed from the scatter within each ensemble. The amplitudes of the mass fluctuations according to 
linear theory, (Tg , at the output times plotted are listed in each panel. The triangles show points for APM 7^ > 20 clusters from Figure 
10 of EDSM. 
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Figure 3. The dependence of the cluster correlation function on the cluster search radius, rc, for the three ensembles at an amplitude 
(Tg = 0.59 (ensembles A and C) and (Tg =1.0 (ensemble B). The cluster catalogues have dc = 35 h-^Mpc. We plot the APM data points, 
as in Figure 2. The points labelled PHG show the correlation function (corrected for projection biases) for the Postman, Huchra &l Geller 
redshift survey of Abell clusters with > (from Figure 10 of EDSM). 
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sidered here, such a large difference cannot be expfained by 
differences in the ciuster sefection. EDSM show that the red- 
shift space correlation function for the PHG sample is highly 
elongated along the line-of-sight, so that fee ~ 1 out to 
line-of-sight separations of ~ 100 /t~^Mpc for projected pair 
separations of < 10 /t~^Mpc. These large anisotropics pro- 
vide strong evidence that the clustering pattern in the R > 
Abell is enhanced by inhomogeneities in the completeness of 
the catalogue. The correction procedure described by EDSM 
removes these biases and brings fee for the Abell sample into 
agreement with the results for the APM sample. 

In summary, for each ensemble we find that cluster cat- 
alogues with similar values of dc have similar two-point cor- 
relation functions. The clustering properties are therefore 
dependent on the shape of power-spectrum of the mass fluc- 
tuations, but not on the way in which clusters are identified. 
This agrees with observational results for APM and Abell 
clusters, but only if the Abell catalogue is corrected for bi- 
ases in the cluster selection. 



3.3 Richness dependence of the correlations 

For each model, and for each of our two cluster selection 
radii, we generated 13 cluster catalogues at each output 
time by varying the lower richness bound so that the mean 
intercluster distance changed from dc ~ 10 /t~^Mpc to 
dc ~ 70 /t~"'Mpc in steps of 5 /t~"'Mpc. We computed the 
two-point correlation function for these catalogues in log- 
arithmic bins using the estimator of equation (3) and we 
computed the radial separation ro at which {cc(''o) = 1 by 
linear interpolation. The dependence of ro against dc pro- 
vides a measure of how the amplitude of the correlation 
functions change with cluster richness. 

In Figure 4, we plot the mean value of ro for each en- 
semble against dc for two values of ag and for the two clus- 
ter selection radii. The error bars on ro show the la errors 
determined from the scatter within each ensemble. As ex- 
pected from the previous subsections, the amplitude of the 
cluster correlation function is insensitive to the amplitude 
of the mass fluctuations and to the cluster detection radius. 
The results for the two F = 0.2 models are virtually in- 
distinguishable, but lie significantly higher than those for 
F = 0.5. For each ensemble, the value of ro becomes ex- 
tremely insensitive to the richness of the cluster sample for 
dc ^ 30 /t~^Mpc. We show in section 4 that the dependence 
of ro against richness is in qualitative agreement with an- 
alytic calculations of the clustering of peaks in a Gaussian 
density field. 

The solid lines in Figure 4 show equation (1), which 
Bahcall & Cen (1992) argue provides a good fit to the ob- 
served cluster correlations. Clearly the strong dependence of 
ro on dc of this relation provides a very poor match to the 
results from our simulations. BC92 have run N-body simu- 
lations of CDM models and seem to find a stronger depen- 
dence of ro on dc than we find from our models for similar 
initial power spectra. In fact, BC92 conclude that equation 
(1) can be reproduced in a CDM model with F Si 0.1. We 
make the following remarks concerning these discrepancies. 
[1] BC92 run a single particle-mesh simulation containing 
250'^ particles for three sets of cosmological parameters in a 
cubical box of side £ = 400 /t~^Mpc. The force resolution of 



their simulations (> 800/t~^kpc) is considerably worse than 
ours, there are some differences in the cluster finding algo- 
rithms, and their mass resolution is better than ours by a 
factor of 6.6. Clearly these differences will lead to differences 
in the structure of the clusters that form in the simulations. 
However, the insensitivity of our results to the cluster selec- 
tion radius described in Section 3.2 suggests that the am- 
plitude of the correlation function is probably not sensitive 
to the resolution of the code. Whereas we have run 10 sim- 
ulations in each ensemble, BC92 run only one model and so 
sample a volume of space smaller than ours by a factor of 
4.2. The errors on their correlation functions will be larger 
than ours and it is plausible that the difference between our 
respective ro-dc relations are caused by statistical fluctua- 
tions. In fact, over the range 10 /t~^Mpc ^ dc ^ 4:0 /t~^Mpc, 
where the errors in ro are small (Figure 4), our results agree 
well with BC92. BC92 use the same phases in each of their 
models (Bahcall, private communication) and this could ex- 
plain why they see a stronger dependence of ro on dc in 
each of their models. From the trends shown in Figure 4, 
a F = 0.1 model should have slightly higher values of ro 
than a F = 0.2 model, but we would expect ro to be nearly 
independent of cluster richness. 

[2] Equation (1) predicts ro w 22 /t"^Mpc for Abell R> 1 
clusters and ro ~ 38 /t~"'Mpc for Abell R > 2 clusters. As 
summarized in Section 1, redshift surveys of Abell clusters 
do indeed give values of ro compatible with these numbers, 
but there is considerable evidence that the Abell catalogue 
contains inhomogeneities that boost the amplitude of the 
clustering pattern. For example, Sutherland (1988) and Ef- 
stathiou et al. (1992) find ro w 14 /t"^Mpc for Abell R>1 
clusters after correction for large scale anisotropics. As we 
have discussed in the introduction, in our view, the empirical 
basis for equation (1) is weak. 

In Figure 5 we show the correlation functions from 
our simulations for two richness cuts matching the space 
densities of 7^ > 20 (dc x 35 h'^Mpc) and 7^ > 35 
(dc X 45 fe-^Mpc) APM clust ers. The observational points 
are from Figure 9 of EDSM. As expected from Figure 4, the 
correlation functions for the models are weakly dependent 
on dc. The APM data are in excellent agreement with the 
two F = 0.2 models, but not with results for F = 0.5. The 
richness dependence of fee for our models is shown in Figure 
7 below for a wider range of dc and compared to calculations 
of the clustering of peaks in a Gaussian density field. 



4 COMPARISON WITH THE STATISTICS OF 
GAUSSIAN DENSITY FIELDS 

In this section, we compare our results with calculations of 
the cluster correlation function using the statistics of Gaus- 
sian density fields. This problem has been tackled by vari- 
ous authors (e.g. Kaiser 1984; BBKS; Bardeen, Bond & Efs- 
tathiou 1987 ). Recently, Mann et al. (1993) have applied an 
analytic technique described by Bond & Couchman (1988) 
to compute the evolution of the correlation function of peaks 
in a Gaussian density field assuming that the peaks move 
according to the Zeldovich (1970) approximation. Holtzman 
and Primack (1993) have used similar techniques to those 
described by Bardeen et al. (1987) to compute the cluster 
correlation function in mixed dark matter models. 
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Figure 4. The dependence of tq for the cluster correlation function on intercluster separation, for two values of (Tg, 0.59 and 1.0, and 
two values of the cluster detection radius, Vc = 0.5 and 1.5 /i~"^Mpc. The solid lines show the relation tq = OAdc. We plot estimates of 
To from the APM and EDCC redshift surveys. 



We use the approximate formula for the peak-peak cor- 
relation function given in BBKS 

ipk.pk = (^^^ + 1) ip.p (5) 

applied to peaks above a sharp threshold i^t- Here (if) is the 
averaged peak height above the threshold Ut and a^(Rs) is 
the variance of the density field smoothed by a Gaussian 
filter of width Rg. To apply equation (5) we need to select 
values for Rs and i^t- Clearly, we should select peaks which 
can be identified with perturbations that will form clusters 
of galaxies; this suggests that we should choose Rs so that 
the mass contained within the Gaussian window, 

M (Rs) = {27r f^^pRl, (6) 

where p is the mean mass density of the universe, is equal 
to the mass of a rich cluster of galaxies. Equation (6) gives 

R,=2.8nl'^( V^"fe"'Mpc. (7) 

The mass within an Abell radius, rA = 1.5 /t~^Mpc (roughly 
the virialized region) of clusters with an abundance equal to 



half that of Abell R > 1 clusters {dc ~ 66 h-^Mpc{EDSMy} 
is estimated to be about 4.2-5.5 x IO^A'^Mq (White et 
al. 1993). The cluster mass functions of our N-body models 
are so steep that if we choose the amplitude of the mass 
fluctuations to approximately match the observed masses 
of clusters with dc = 63 /t~^Mpc, then we find that the 
masses of clusters with the APM 7?. > 20 abundance (dc = 
35 /t"^Mpc) are f» 2 X lO"/t"^M0 (see Figure 1 of White 
et al. al. (1993)). Thus, equation (6) gives about the correct 
masses for clusters over a wide range of abundances if we 
set Rs K. 5QJ^^ fe-^Mpc. 

Having chosen Rs, we can fix the threshold Pt by match- 
ing the observed number density of clusters to the number 
density of peaks obtained by integrating the differential peak 
density distribution (BBKS equation 4.3) 

M,u{u)du= I 3 e"^G(7,7^), (8) 

where G is a function given in Appendix A of BBKS and 
the parameters 7 and i?, are given by moments of the power 
spectrum. 

This method of identifying peaks with real clusters is 
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Figure 5. The correlation function for clusters identified with Vc = 0.5 h "^Mpc with same mean separation as APM clusters with 
richness K > 20 (dc = 35 /i-^Mpc) and K > 35 (dc = 45 /i-^Mpc). 



clearly very rough, so it is useful to compare the results with 
another method of choosing Rs and i^t- Instead of fixing Rs 
we choose Rs and i^t to reproduce the desired abundance of 
clusters and to satisfy 

i^MRs) = Sc, (9) 

where 8c, the critical density contrast for forming virialized 
objects by the present day, is set equal to 1.69 independent 
of Qo (see White et al., 1993). This scheme for choosing Rs 
and Pt has been used by BBKS and Mann et al. (1993). 

Figure 6 shows calculations of the correlation length as 
a function of dc- We have adopted erg = 0.59 for the models 
A and C, and erg = 1 for model B so that the figure can 
be compared with the results from the N-body simulations 
plotted in Figure 4. For the Q = 1 models, the two methods 
for choosing Rs and Pt give similar results at large values of 
dc- However, for ensemble B (Qo = 0.2), equations (8) and 
(9) give Rs closer to 5 /t"^Mpc at dc > 40 /t"^Mpc than the 
value Rs = 2.5 /t"^Mpc that we used to compute the dot- 
ted line. Comparing Figures 4 and 6 we see that the simple 
theory described here consistently overestimates the ampli- 
tude of the correlation function measured in the simulations. 
Mann et al. found a similar disagreement in comparing their 
theoretical predictions with the N-body results of BC92. 

The discrepancies between Figure 6 and Figure 4 are 
worse for the higher value of erg , and are of the same order 
as the correction for the dynamical evolution of the mass- 
density fluctuations in equation (5) (approximated by the 
addition of to p/cr). It might be thought that a more so- 
phisticated treatment of the dynamical evolution of the cor- 
relations could improve the agreement. However, the com- 
putations of Mann et al. (1993), who use the Zeldovich ap- 
proximation to model the motions of peaks, give a worse 
match to our results than the simple theory outlined above. 



The dynamical contribution to the clustering is significant 
even in the simple analytic examples plotted in Figure 6 
(see also Figure 7), i.e. iz/cr ~ 1 and it is unclear how well 
the Zeldovich approximation describes the non-linear clus- 
tering in the numerical simulations. Furthermore, the peak 
prescription provides only a very crude guide to the loca- 
tions and masses of clusters (see e.g. Bond et al. 1991). It is 
not particularly surprising therefore that the more elaborate 
calculations of Mann et al do not seem to provide a better 
description of the N-body results than the simple model de- 
scribed above. Evidently, a more precise description of the 
non-linear evolution of the mass fluctuations, possibly in- 
cluding cluster merging, is required to describe the N-body 
results. The qualitative trends shown in Figure 6 do match 
the N-body results; the relative amplitudes of fee for the 
three models, and the levelling-off of ro with increasing dc, 
are in rough agreement with the results described in Section 
3. 

In Figure 7, we plot the correlation functions for three 
cluster abundances at a single output time for each of the 
three ensembles. The dotted lines show the predictions of 
equation (5) with Rs and Pt fixed by equations (8) and (9). 
In each case the dotted lines overestimate the amplitude of 
the N-body curves. The shapes of the correlation functions 
for r = 0.2 are quite well described by the linear theory 
shape of {p,p, whereas the N-body results for F = 0.5 are 
noticeably steeper than the linear correlation function on 
scales < 10 /t"^Mpc. 

In summary, simple models based on the statistics of 
high peaks can provide a rough match to the clustering of 
rich clusters measured in our simulations. However, the am- 
plitude of fee is consistently overestimated by the high peak 
model. 
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Figure 6. The dependence of the correlation length tq against intercluster distance dc computed from equation (5) for the three ensembles 
A-C. The solid lines show computations in which Rg and vt satisfy equation (9) and the dotted lines show results for a fixed value of Rg 
(listed in /i~^Mpc in each panel). The dashed line shows the Bahcall Cen relation (equation 1). 




Figure 7. Cluster correlation functions for the three ensembles at (Tg = 0.59 (ensembles A and C) and (Tg = 1.0 for three values of dc-, 
dc = 15,35 and50 /i-^Mpc. The errors on the data points show one standard deviation of the mean. The dotted lines show the predictions 
of equation (5) in which Rs and ut satisfy equation (9). Clusters in the N-body simulations were identified with Tc = 0.5 h~^Mpc. 



10 R. A. C. Croft and G. Efstathiou 



5 SIMULATED APM REDSHIFT SURVEYS 

In this section, we use our simulations to generate cluster 
catalogues with the same selection criteria and volume as 
the APM cluster redshift survey of Dalton et al. (1992). We 
generate the simulated catalogues in the following manner: 

We select rich cluster catalogues from the simulations 
normalised to the mean space density of the APM TZ > 20 
clusters and we replicate our computational box out to a 
redshift of 50000km s~^ ( 5/3 boxlengths) using the periodic 
boundary conditions. After projecting the cluster positions 
onto the surface of a sphere, we discard all those outside the 
APM survey area. We select clusters at random using the 
selection function for APM clusters (see EDSM Figure 1) 
resulting in simulated catalogues with statistically similar 
number-redshift distributions and total number of clusters 
as the real survey. The two point cluster correlation function 
in redshift space for the simulated catalogues was computed 
by cross-correlating with a random catalogue and using the 
estimator of equation (1) from DEMS. 

Our results are shown in Figure 8. The lines in each 
panel show fee averaged over 10 simulations of the APM sur- 
vey, and the error bars show the la variation of the mean. 
As above, the F = 0.2 models provide an excellent match to 
the observations. The cluster random peculiar velocities in 
the simulations average ~ 600km s~^, and so redshift space 
distortions have only a small effect on the cluster correla- 
tion functions on separations greater than > 5 /t~^Mpc (c/. 
Figures 2 and 6). At smaller separations, redshift-space dis- 
tortions lead to a reduced amplitude for fee, which agrees 
better with the observations than the real-space correlation 
functions plotted in Figure 2. The APM cluster redshifts 
have errors of ~ 500km s~^ which have not been modelled 
in the simulated redshift surveys. The amplitudes of the cor- 
relation functions would have been reduced slightly further 
on small scales had we included redshift errors in the simu- 
lated catalogues. 



Table 2. Errors on 
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0.55 
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0.89 
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0.60 



The error bars on the APM data points are 'Poisson' 
errors computed from the formula 



(10) 



+ 

where Ncc is the number of distinct cluster pairs in the bin 
centred at separation s in redshift space. Our simulations 
show that equation (10) usually underestimates the errors 
on fee. In Table 2 we give the mean error on fee for a single 
simulation estimated from equation (10) ((5{cc)p) divided 



by the scatter measured from the 10 simulations ((t({cc)) 
as a function of radius for the three ensembles. Equation 
(10) underestimates the errors in ensemble A by a factor of 
about 1.1, on scales s > 10 /t~^Mpc. For ensembles B and C, 
equation (10) underestimates the errors by a factor of about 
1.5-1.7 for s > 10 /t"^Mpc.This is in rough agreement with 
the error analyses of Ling, Frenk & Barrow (1986) and Mo, 
Jing & Borner (1992) based on bootstrap resampling applied 
to real samples. Although the errors are underestimated by 
equation (10), none of our simulated APM surveys with F = 
0.5 match the high amplitude of fee for the real survey at 
s > 8 /t"^Mpc. Thus the standard CDM model gives a poor 
description of the APM survey. 



6 CONCLUSIONS 

We have described a series of TV-body simulations of CDM- 
like models designed to study the clustering of rich clusters. 
We find the following results: 

[1] The cluster correlation function is insensitive to the am- 
plitude of the mass fluctuations. 

[2] The cluster correlation functions for catalogues con- 
structed with two values of the cluster radius Tc, but with 
the same intercluster separation, are almost identical. The 
correlation functions are thus insensitive to the cluster selec- 
tion, for values of Tc comparable to those used in real cluster 
samples. 

[3] In each set of simulations the amplitude of fee is found 
to depend weakly on cluster richness. For clusters with 
dc > 30 /t"^Mpc, we find correlation lengths of ro ~ 9- 
10 h-^Mpc for the standard CDM model, (F = Q/t = 0.5) 
and ro ~ 14 -15 /t~^Mpc for low density CDM models with 
F = 0.2. 

[4] The cluster correlation functions in the low density F = 
0.2 models are almost identical to those in the Q = 1, F = 0.2 
models. The cluster correlation function in spatially flat 
models is therefore insensitive to the value of the cosmolog- 
ical density parameter but strongly dependent on the shape 
of the initial fluctuation spectrum. 

[5] The high-peak model for fee described in Section 4 over- 
estimates the amplitudes of the correlation functions mea- 
sured in our models, but qualitatively reproduces the ob- 
served trends with richness and with the shape of the power 
spectrum. 

Comparing the results from our simulations with obser- 
vations of the clustering of rich clusters we conclude: 
[6] The standard CDM model provides a poor match to the 
clustering measured in the APM and other cluster redshift 
surveys. The APM survey is well described by a low density 
CDM model with F w 0.2. 

[7] The weak dependence of ro on dc seen in our models does 
not agree with results for the clustering of rich Abell clusters. 
In particular, the very high amplitude ro = 45 ± 5 /t~^Mpc 
measured by Peacock and West (1992) for Abell R > 2 clus- 
ters is difficult to reconcile with the CDM-like models con- 
sidered here. A similar conclusion has been reached by Mann 
et al. (1993). 

The observed clustering of galaxies agrees well with the 
power-spectrum of a low density, F Si 0.2, CDM model (see 
Efstathiou, Sutherland & Maddox 1990, Vogeley et al. 1992, 
Loveday et al. 1992). Our results suggest that the clustering 
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Figure 8. The correlation function for simulated catalogues of clusters compared to the APM data points. We plot the mean of 10 
simulated catalogues in each panel and show the 1(T scatter of the mean. The errors on the APM points are the Poissonian errors 
computed from equation (10). 



of rich clusters selected from digitised photographic plates 
is also consistent with such a model. However, there is po- 
tentially a discrepancy with the clustering of richer clus- 
ters selected from the Abell catalogue. If, as Bahcall and 
West (1992) and other authors have argued, rich clusters 
follow a relation of the form ro = OAdc, then we would 
probably have to abandon theoretical models for the origin 
of structure based on Gaussian initial conditions. However, 
the clustering of rich Abell clusters may be high because of 
incompletness in the Abell catalogue, and as Peacock and 
West note, their quoted error on ro for R > 2 could be an 
underestimate. Clearly it would be interesting to complete 
a redshift survey of several hundred rich clusters from one 
of the automated cluster catalogues. According to our nu- 
merical simulations, and the observational results already 
obtained for the APM and EDCC clusters, we would not 
expect ro to rise much above x 16 /t~^Mpc. 
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